Let \(\mathcal{F}\) a parametric statistical model for the data \(y\), where \(f(y;\theta)\) is the density or probability function and \(\theta \in \Theta\) is a p-dimensional parameter. By considering \(f(y;\theta)\) as function of \(\theta\) with \(y\) fixed to the observed value, then the likelihood function of \(\theta\) based on \(y\), \(L: \Theta \rightarrow \mathbb{R}^{+}\), is \[L(\theta) = L(\theta; y) \propto f_Y(y; \theta) \,\,, \] where the proportionality is related to a constant \(c(y)>0\) that does not depend on \(\theta\).

On the basis of the data \(y\), \(\theta \in \Theta\) is more likely than \(\theta' \in \Theta\) as an index of the data generating model if \(L(\theta) > L(\theta')\)

Let \(y=(y_{1},\ldots, y_{n})\) a sample of i.i.d. values from a Bernoulli distribution, \(Y \sim \mbox{Be}(p)\). Then the likelihood function is

\[\mathcal{L}(p) = \prod_{i=1}^{n} p^{y_i} (1-p)^{1-y_i} \] and the log-likelihood function takes the form \[\ell(p) = \sum_{i=1}^{n} y_i \log (p) + (1-y_i) \log (1-p)\] The maximum likelihood estimate, \(\hat p\), is the sample proportion \(\hat p =\frac{1}{n}\sum_{i=1}^{n}y_i\). To derive it we obtain the score function \[\ell_{\star}(p)=U(p)=\frac{\partial \ell(p)}{\partial p}= \sum_{i=1}^{n} \frac{y_i}{p} - \frac{1-y_i}{1-p} = \frac{\sum_{i=1}^{n}y_i}{p} - \frac{n- \sum_{i=1}^{n}y_i}{1-p}\] By equating at zero the score function, we obtain the maximum likelihood estimate

\[U(p) = 0 \implies \hat p = \frac{1}{n}\sum_{i=1}^{n} y_i \]

Suppose generating a random sample \(y\) of size \(n=100\) from Bernoulli distribution with parameter \(p=0.6\). We want make inference on \(p\). Thus,

We want write a function taking in argument the parameter and the data and returns the log-likelihood. We can do it in two ways:

• Leveraging the dbinom() function

llik_bin <- function(theta, data){
  sum(dbinom(x = data, size = 1, prob = theta, log = TRUE))
}

• Writing the log-likelihood function manually

llik_bin2 <- function(theta, data){
  sum(data) * log(theta) + (length(data) - sum(data)) * log(1-theta)
}

We can check if both the functions return the same value. Thus, for instance by fixing \(p=0.5\)

llik_bin(0.5,y)

[1] -69.31472

llik_bin2(0.5,y)

[1] -69.31472

The maximum likelihood estimate is in closed form (later we will see how carrying out the ML using numerical optimisation, which is particularly useful when the ML estimate cannot be obtained analitically), so

MLEp  <- mean(y)
MLEp

[1] 0.63

llik_bin(MLEp, y)

[1] -65.89557

pgrid <- seq(0.01, 0.99, by = 0.01)

pgrid[which(llik_bin2(pgrid, y) == max(llik_bin2(pgrid, y)))]

[1] 0.63

Note that above we used the function built manually. This because that function is vectorized, while the first one should be vectorized before using it. Indeed

llik_bin(c(0.5, 0.75, 0.99), y) # Wrong

[1] -98.22045

llik_bin2(c(0.5, 0.75, 0.99), y) # Correct

[1]  -69.31472  -69.41686 -171.02447

# The results above is due to this operation
sum(dbinom(x = y[seq(1,100,3)], size = 1, prob = 0.5, log = TRUE)) + 
sum(dbinom(x = y[seq(2,100,3)], size = 1, prob = 0.75, log = TRUE)) + 
sum(dbinom(x = y[seq(3,100,3)], size = 1, prob = 0.99, log = TRUE))

[1] -98.22045

# However, we can vectorize the funtion llik_bin by means of 
llik_bin_v <- Vectorize(llik_bin, 'theta')
llik_bin_v(c(0.5, 0.75, 0.99), y) # Correct

[1]  -69.31472  -69.41686 -171.02447

Then we can plot it, including two vertical lines denoting where the true parameter value and the maximum likelihood estimate are located

par(mfrow = c(1, 2))
curve(llik_bin2(theta=x, data = y), 0, 1, 
      xlab= expression(p), ylab = expression(l(p)),
      main = "Log-likelihood function")
abline(v = p, col = "red")
abline(v = MLEp, col = "blue")
legend("topleft", legend = c(expression(p), expression(hat(p))),
       col = c("red", "blue"), lty = c(1, 1))
curve(llik_bin_v(theta=x, data = y), 0, 1, 
      xlab= expression(p), ylab = expression(l(p)),
      main = "Log-likelihood function")
abline(v = p, col = "red")
abline(v = MLEp, col = "blue")
legend("topleft", legend = c(expression(p), expression(hat(p))),
       col = c("red", "blue"), lty = c(1,1 ))

We know that \(\hat \theta \stackrel{\cdot} \sim \mathcal{N}(\theta,\mathcal{I}^{-1}(\theta) )\), where \(\mathcal{I}(\theta)\) is the expected information matrix, that is \(\mathcal{I}(\theta)=E(J(\theta; Y))\), with \(J(\theta; Y)\) the observed information matrix, that is the negative hessian. In our one-dimensional parameter example \[J(p) = -\frac{\partial^2 \ell(p)}{\partial p^2} = \bigg[\frac{\sum_{i=1}^{n}y_i}{p^2} + \frac{(n-\sum_{i=1}^{n}y_i)}{(1-p)^2}\bigg]\] Thus, since \(\sum_{i=1}^{n} {Y_i}\sim Bin(n,p)\) we have \[\mathcal{I}(p)=E[J(p; Y)]= \bigg[\frac{\sum_{i=1}^{n}E[Y_i]}{p^2} + \frac{(n-\sum_{i=1}^{n}E[Y_i])}{(1-p)^2}\bigg] = \frac{n}{ p (1-p)}\] Then \(\hat p \stackrel{\cdot}\sim \mathcal{N}(p, \frac{p(1-p)}{n})\), where we can replace \(\mathcal{I}^{-1}(p)\) with the estimate \(i^{-1}(\hat p)=j^{-1}(\hat p)=\hat p (1- \hat p)/n\)

Let us consider the logit function \(\psi(p) = {\rm logit}(p)=\log\bigg(\frac{p}{1-p}\bigg)\), which gives the log-odds. We can use it to reparametrise the model, so the parametric space is unbounded. At such point, the parameter is expressed in the logit scale, and we need to re-express them in the original scale, that is \(p(\psi) = \exp(\psi)/(1+\exp(\psi))\). We can leverage the property of invariance under reparametrisations of the maximum likelihood estimator and so the maximum likelihood estimate is simply \(\hat \psi = \log\bigg(\frac{\hat p}{1- \hat p}\bigg)\). Then, we do not need to write down the log-likelihood under the reparametrisation obtain the MLE of \(\psi\) from it.

However, we can easily visualize the log-likelihood function under the reparametrisation as

The asymptotic distribution of the ML estimator of \(\psi\) is given by \(\hat \psi \stackrel{\cdot} \sim \mathcal{N}(\psi, V(\hat \psi))\)

By using the delta method, we can easily obtain

\[ V(\hat \psi) = V(\hat p) \bigg(\frac{d}{d p} \psi (p) \bigg)^2 = V(\hat p) \bigg(\frac{d}{d p} \log\frac{p}{1-p}\bigg)^2 = \frac{ p (1 - p)}{n} \frac{1}{p^2(1-p)^2} = \frac{1}{np(1-p)}\] and a consistent estimate of such variance is \(\hat V(\hat \psi) = \frac{1}{n \hat p (1 - \hat p)}\)

In this example, we computed analytically some quantities and we obtain numerical values just by plugging into the formulas the input. However, remind the MLE for a parameter may not have an analytic solution. Many times we do not have a closed form for MLE estimates and we may need numerical optimisation. R provides various functions for performing numerical methods:

• nlm(): minimizes a function using a Newton-type algorithm. It needs a starting value and does not allow constraints on the parameters. It is usually fast and reliable. It returns the ML estimate \(\hat{\theta}\) (estimate), the value of the likelihood \(-l(\hat{\theta})\) (minimum) and the hessian (hessian), if hessian = TRUE.

• optim(): minimizes a function using Nelder-Mead, quasi-Newton and conjugate gradients algorithms. It includes an option for box-constrained optimization, and it requires a starting value. It returns the ML estimate \(\hat{\theta}\) (par) and the value of the likelihood \(-l(\hat{\theta})\) (value) and the hessian (hessian), if hessian = TRUE. You also can maximize the function by using fnscale = -1 in control argument.

• nlminb(): often is more stable, robust and reliable than optim (in particular with “nasty” functions). It performs only minimization and does not yield numerical derivatives as output. It returns the ML estimate \(\hat{\theta}\) (par) and the value of the likelihood \(-l(\hat{\theta})\) (objective).

• optimize(): by using a combination of golden section search and successive parabolic interpolation, searches in an interval for a minimum or a maximum (if maximum=TRUE) of a function. It returns the ML estimate \(\hat{\theta}\) (minimum) and the value of the likelihood \(l(\hat{\theta})\) (objective). Drawback: suited only for one-dimensional parameter.

Here, we are simply using the Binomial example to show how (some of) these functions work.

Note: In the following, we consider and implement the negative log-likehood function, since some numerical optimisers that we will see are only able to perform minimisation of functions.

In addition, we implement the score and the observed information (used for comparison with the optimiser output)

To check if the score and the observed information are correctly implemented, we can compute the numerical gradient and hessian of the log-likelihood function, as implemented in the {numDeriv} R package.

Let’s start with the first optimisernlm(considering as starting value 0.5, which is quite close to the MLE; however, see the warnings printed on the console)

The previous optimiser does not allow constraint on the parameters; however, there are some replacements in the algorithm that allow you to reach the result. If you want optimise a function which is constrained it should be better to consider optimiser allowing for box constrained optimisation, as for instance optim.

• The optim() function provides a lot of numerical methods (Nelder-Mead, quasi-Newton, conjugate-gradient methods and simulated annealings). As a drawback, the user has to set up carefully the initial guesses and the method, since the final solution may be quite sensitive to these choices…To compute the observed information, the function optimHess() computes numerical derivatives of generic functions.

• As an alternative, we could consider again the reparametrisation introduced above. Generally, for numerical purposes, it is convenient to reparameterize the model in such a way that the new parameter space is unbounded. Then, we explore this by using nlm()

Regression is trying to describe a dependent variable \(\mathbf{y}\) through some statistically significant predictors \(\mathbf{x}\), linking them through a function \(f(\cdot)\) (also referred as systematic component): \[\mathbf{y} = f(\mathbf{x})+\boldsymbol{\varepsilon},\] where \(\boldsymbol{\varepsilon}\) is a measure of noise implicit in modelling (also stochastic component).

Checking the assumptions of a regression model is not just a theoretical matter, but it is something related to scientific reproducibility, scientific transparency and visualization power. Although it is impossible to write down an exhaustive list of ‘what to do for checking a model’, we put here just some initial steps for visualising and inspecting your variables, as well as to inspect possible correlations between them and the inclusion of further predictors. The suggested steps refer to the (normal) classical linear model, \(f(\mathbf{x}) = f(\mathbf{x}; \boldsymbol{\beta}) = \mathbf{X} \boldsymbol{\beta}\), that is

\[y_i=\beta_1x_{i1} + \beta_2x_{i2} + \ldots + \beta_px_{ip} +\varepsilon_i, \quad \varepsilon_i \stackrel{iid}\sim \mathcal{N}(0,1)\] where we usually assume \(x_{i1}=1\), but they contain principles that can be extended to statistical modelling in general.

• Examine the distribution of each of the explanatory variables, and of the dependent variable. Thus, skewness and presence of outliers can be detected. When the distribution is skew, consider transformations inducing symmetry.
• Examine the scatterplot involving all the explanatory variables. If some pairwise plot shows non-linearity, consider some possible transformations.
• Examine the range of the variables.
• Examine the degree of correlation between the explanatory variables.

• Plot residuals against fitted values and check for patterns in the residual. Also, check for homoscedasticity.
• Examine the Cook’s distance for possible outliers.

We consider the dataset nihills available in the {DAAG} package, containing data from the 2007 calendar for the Northern Ireland Mountain Running Association. The number of observations (races) is \(n = 23\), and for each of them we register the following variables:

• Distances in miles: dist.
• Amount of climbs expressed in feet: climb.
• Record time in hours for males: time.
• Record time in hours for females: timef.

Objective: We want to regress the times through some predictors.

In the following, we focus on the time results for males. Then, we load the data.

We explore the data structure, the first six and the last six rows.

A summary report of the variables allows exploring some characteristics of the variables.

In a first step, we consider only the predictor dist for regressing the dependent variable time. Let us have a quick look to the data by means of a scatterplot:

Apart one, all the data times are lower than 2 hours. As an initial attempt, let’s fit a simple linear model (and see a basic model output): \[\texttt{time}_i =\beta_1+ \beta_2\texttt{dist}_i+\varepsilon_i, \ \ i=1,\ldots,n, \quad \varepsilon_i\stackrel{iid}\sim \mathcal{N}(0,\sigma^2).\]

A more detailed summary of the fitting procedure can be obtained by using the summary() function

• Call: The fitted model.
• Residuals: Some residual summaries, mainly quantiles (asymmetry?).
• Coefficients:
  • Estimates: coefficient estimates, \(\hat \beta_j, j = 1, \ldots, p\).
  • Std. Error: estimated standard deviation of the estimators \(\hat \beta_j\), that is \(\sqrt{ \widehat {V (\hat \beta_j)}}\).
  • t-value: the (observed) values of the test statistic for testing the nullity of the coefficients, \(H_0: \beta_j=0\) against \(H_1: \beta_j\neq0\), that is \(t_j = \hat \beta_j / \sqrt{ \widehat {V (\hat \beta_j)}}\).
  • \(Pr(>|t|)\): the corresponding p-values, \(2P(T_{n-p} > |t_j|)\).
  • Some symbols denoting where the p-values fall, *** if \(0 < p \leq 0.001\), ** if \(0.001 < p \leq 0.01\), * if \(0.01 < p\leq 0.05\), . if \(0.05 < p \leq 0.1\), otherwise nothing.

• Residual standard error: the square root of the unbiased sample variance of the residuals, that is \(\sqrt{s^2}\), with the corresponding degrees of freedom \(n-p\).
• Multiple R-squared: the \(R^2\) coefficient.
• Adjusted R-squared: the adjusted \(R^2\) coefficient (given by \(R^2 -\frac{p-1}{n-p}(1-R^2)\)).
• F-statistic: the observed value for the test statistic comparing the full model with the null model (that is the model including only the intercept), with the corresponding degrees of freedom and the related p-values. Namely, for testing \(H_0: \beta_2= \ldots = \beta_p=0\) against the alternative that at least one is different from zero.

Then, we can visualise some aspects of the fitted model. The first step, and likely the most important one, is to explore the residuals plots

• (Raw) residuals: \(e_i = y_i - \hat y_i\)
• Standardised residuals (internal studentized residuals): \(r_i = e_i/(s \sqrt{1-p_{ii}})\), where \(p_{ii}\) is the \(i\)-th diagonal element of the projection matrix (\(P = X (X^\top X)^{-1} X^\top\))

Plots (graphical assessment of certain linear model assumptions):

• Residuals vs fitted: \((\hat y_i, e_i)\) (linearity and homoscedasticity assumptions).
• Q-Q residuals: \((\Phi^{-1}((i-0.5)/n), r_{(i)})\) (normality assumption).
• Scale-Location: \((\hat y_i, \sqrt{|r_i|}\) (homoscedasticity assumption).
• Residuals vs Leverage: \((p_{ii}, r_i)\), with the inclusion of contour level of \(q^2(D,h)=D p (1-h)/h\), that is the lines \(q = \pm \sqrt{D p (1-h)/h}\), where \(D\) is the Cook’s distance.

Then, we can produce the scatterplot including the regression line

Of course, the fit seems good enough. However, the relationship between time and dist seems not to be purely linear…and residual plots seem to suggest a lack of fit for a few points.

Here, we obtain the quantities of the output of the summary() function applied to an lm object by hand. In the following, the computations are done using matrix algebra (albeit we are working with a simple linear model). In this way, the following lines of code can be used to (eventually) carry out some checks using the multiple linear model.